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Abstract 

We exploit the analogy between dynamics of inertial particle pair separation in a random-in-time 
flow and tlie Anderson model of a quantum particle on the line in a spatially random real-valued 
potential. Thereby we get an exact formula for the Lyapunov exponent of pair separation in a special 
case, and we are able to generalize the class of solvable models slightly, for potentials that are real up 
to a global complex multiplier. A further important result for inertial particle behavior, supported 
by analytical computations in some cases and by numerics more generally, is that of the decay of 
the Lyapunov exponent with large Stokes number (quotient of particle relaxation and flow turn-over 
time-scales) as St~'^^^. 
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1 Introduction 



We propose in this paper to study evolution of inertial particle pair separation in a smooth random 
flow. The dynamics of pair separation, through a very simple change of variables, is equivalent to 
the Anderson localization problem jHj of a quantum particle in one-dimensional space with random 
potential, with the twist that the potential is not necessarily a real-valued function. This analogy 
permits us to calculate explicitly, in certain situations, the Lyapunov exponent of inertial particle 
pair separation. In our slightly generalized Anderson problem with complex potential, we find some 
additional cases when the Lyapunov exponent may be computed analytically. The article treats 
the topic from the viewpoint of inertial particles, but some results could also be interesting for the 
Anderson localization problem. 

Understanding the behavior of inertial particles, and more specifically aerosols (tiny liquid or solid 
particles in a gas), has several important applications. For example rain formation inside clouds 
depends on collisions of water droplets jH], and a better understanding of this process is needed 
for reliable weather forecasts. We may also mention internal combustion engines (oil droplets, coal 
powder) where knowledge of particle behavior helps to optimize design. Since we will mostly study 
here the two dimensional situation, let us also mention the problem of particles floating on the free 
surface of a liquid, such as the surface of the ocean 

The Lyapunov exponent of particle separation is one of the most fundamental quantities describing 
relative motion of close-by particles. Finer details require knowledge of the distribution of finite-time 
Lyapunov exponents, which would permit to predict such things as collision rate of particles [2], to 
give an example. This paper is a first step in this direction. 

Let us outline the plan of the paper. First we give a presentation of the smooth Kraichnan flow 
in which our particles shall evolve. We then describe dynamics of inertial particles and in particular 
of pair separation, which we reduce to a two-dimensional problem. After that we treat the special 
case of genuinely 2D flows. Later we exhibit some additional solvable cases which generalize that 
of the classical Anderson model. We then come to results of numerical simulations that sustain our 
findings and finally discuss inexactness of Piterbarg's formula in ^01 as compared with numerics. 

This paper also gives a review of the topic from the point of view of inertial particles, so that 
results already known are rederived within this framework and with our notations, thus serving 
possibly as reference for further investigation in this direction. Effort was made to give detailed and 
clear calculations. 

The author would like to thank Krzysztof Gaw^dzki for sharing his ideas on the problem and 
pointing out the relationship with the Anderson model. Lot of credit goes to Alexander Fouxon, who 
collaborated in this work, details of which are to appear in a common paper j^; this collaboration 
was made possible thanks to the European Network: "Fluid mechanical stirring and mixing: the 
Lagrangian approach". The author is grateful to Angelo Vulpiani who offered the postdoc position 
and grant (funded by COFIN 2003: "Statistical description of complex systems and systems with 
many electrons") while this work was prepared at the TNT group at the Physics Department of the 
Universita di Roma "La Sapienza". 

2 Smooth Kraichnan velocity field 
2.1 The strain matrix 

We consider a smooth Kraichnan velocity field. At small scale it may be linearized so that velocity 
increments behave like 

U{R + r,t) - U{R,t) = {f.V^)U{R,t) = rjdjUi{R,t)ei 

(cj are unit vectors of the chosen basis. Also note the convention that we use uppercase R, U for 
absolute position and speed while later we will use lowercase r, v for the relative position and speed 
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of two particles.) It is convenient to introduce the velocity gradient matrix, also called strain matrix: 

a,,{R,t) = djU^{R,t) (1) 
in order to write U{R + r,t) - U{R, t) = a{R, t)r. 

2.2 Strain matrix correlation tensor 

In the Kraichnan model the velocity field [/ is a Gaussian, white-in-time vector field, also called 
Brownian vector field. We suppose a{R,t) to be a mean zero process (for all R,t), so that it is 
completely determined by its second moment {a{R,t) (g) a{R' ,t')), which is proportional to 6{t — t'), 
because of the white-in-time hypothesis. 

What will be important for us in tracking particle pair separation is the statistics of the strain 
matrix along the trajectory of the reference particle, as can be seen from (fT6|l . Suppose that the 
reference particle follows the trajectory R{t), which we shall always suppose continuous. Then define 

ait)=a{R{t),t) (2) 

Since the matrix field a is centered Gaussian delta-correlated in time, the same will hold for (T{t). 
Thus it is enough to know the equal-time two-point function of a{t) to completely characterize it. In 
fact we have 

{a{t) a{t')) = {a{R{t),t) a{R{t'),t')) 

= {a{R{t),t)(g)aiRit),t')) (3) 

where the first equality follows from the definition ^ of cr(t), and the second equality holds because 
the correlator on its left hand side being proportional to 6{t — t'), one may replace R{t') by R(t), 
since we suppose R{t) continuous in t and we suppose the two-point correlator of the field a to be 
continuous in the spatial variable. 

Let us now introduce the equal-time equal-position correlation tensor Cikji of the field cr, which 
may be defined through the equality 

{a^R, t)ajiiR, t')) = Sit - t') Cik,ji (4) 

Since we suppose the statistics of the field a homogeneous in space and time, Cikji is a constant 
independent of position and time. In Sect. 12.61 we point out that for the velocity field U it is enough 
to have spatially homogeneous increments, without U itself being spatially homogeneous, for a to be 
homogeneous in the position variable. 
Combining JSJ and (@J we have 

Mt)a,iit')) = 6it-t')Cik,ji (5) 

In the setup of the Kraichnan model, where a{t) is Gaussian white noise in time, it is often better 
to introduce and use 

dS = a{t)dt (6) 

and then we can talk about the covariance (better called covariance process in the mathematical 
literature) of dS and write the more correct 

{dSik,dSji) = Cikjidt 

The covariance of cr(t) may be derived from the velocity two-point function {Ui{R,t)Uj{R' ,t')) as 

{aik{R,t)aji{R,t')) = {{dkUi{R,t)){diUj{R,t'))) (7) 
= dR,dn,mR,t)U,iR',t'))\^,^^ (8) 
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2.3 Translation invariance 



Supposing spatial homogeneity of the velocity statistics, which may be expressed as 

{Ui{R,t)UjiR',t')) = 6it - t')D,j{R-R') 
for some tensor-valued function D, permits to further write ((HI as 

{aik{R,t)aji{R,t')) = -6{t - t'){dkdiD,,){0) 
which in particular, upon substitution into Q, leads to the relationship 

Cik,ji = -{dkdiDij){6) (9) 

We thus arrive at the important conclusion that if the statistics of the velocity field U is spatially 
homogeneous, then Cikji is symmetric under exchange of k with /, or equivalently of i with j, since 
the symmetry under simultaneous exchange of both pairs is granted by the definition Q. 

Note however that one can imagine flows where, at least in some bounded region of space, statistics 
of a is homogeneous but that of the velocity field U is not. For example a uniformly expanding flow 
is like that. To give examples of such a flow, one can imagine either the surface flow of some liquid 
with adequate up- and down- welling in some region, or the surface of a bubble which grows and 
shrinks in time (but that is a curved surface, and we do not pretend to treat that situation in this 
paper), or perhaps a flat film whose boundary is attached to a growing and shrinking ring. 

2.4 Isotropy 

Coming back to the case when the velocity field is spatially homogeneous, if furthermore we sup- 
pose it isotropic and not breaking parity symmetry (defined here as orthogonal transformations of 
determinant —1), we have the second-order development for D (the velocity field is smooth!) 

Ai(r) « DqS.j -Di[{d+1- 2p)5,jr^ + 2{pd - ly.r^] + o{r^) 

where p is the compressibility degree 0] of the flow and Z?i is a dimensional constant of dimension 
time"^. Hence, from © 

Cik,ji = -{dkdiDij){0) = 2Di[{d + 1 - 2p)6ij5ki + {pd - I)i5ik6ji + 5udjk)] 

2.5 Diagonalization, positivity 

As a side note, it is interesting to see that Cikji is indeed a positive matrix, and if we wanted to 
generate a from independent white noises, it would be useful also to be able to diagonalize Cikji- 
We start with the latter, and criteria for positivity will be obvious then. 

Consider Cikji as the matrix of a linear transformation acting on d x d matrices by 

(^ji '-^ (^ik = Cik,ji<^ji (10) 

We identify three terms in Cikji, namely Sij6ki, ^ik^ji and SuSjk- The first one is clearly the identity, 
which we shall conveniently denote Id. The second acts as M i— > (trM)!^, (where Id is the dx d 
identity matrix, that is diag^(l, . . . , 1)) and we shall call it Tr. The third acts as transposition and 
we will call it Tp. 

Prom this it is clear that the three operators Id, Tr and Tp commute mutually and thus may be 
diagonalized in a common basis. It is also straightforward to find this basis: use as eigenspaces the 
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space 'asym' of anti-symmetric matrices, the space 'symO' of symmetric traceless matrices, and 'seal' 
of scalar matrices. We have the following eigenvalue structure: 





Id 


Tr 


Tp 


Cik,jl 


asym 


1 





-1 


2Di{d + 2){l-p) 


symO 


1 





1 


2Di[{d-2){p + l) + 2\ 


seal 


1 


d 


1 


2Dip{d + 2){d-l) 



For Cik,ji to be positive on each eigenspace we get thus, supposing d > 2 and Di > 0, the 
necessary and sufficient condition < p < 1. 

Since the correlation matrix Cikji is symmetric, it can be diagonalized in a basis which is or- 
thonormal with respect to the scalar product induced by the identity matrix, which here means 
^ij^ki (see that this tensor is equal to 1 exactly when the ordered pairs (i, k) and {j, I) coincide, and 
equal to otherwise), and this scalar product is simply the sum of products of corresponding matrix 
elements, formally given by {cr,a') = (Jik(T[;^- In order to generate a (vector-valued) Gaussian random 
variable with covariance matrix Cikji, it is enough to multiply each element (an eigenvector) of a 
given orthonormal diagonalizing basis of Cikji by the square-root of the corresponding eigenvalue 
and by an independent (for each basis element) standard normal random variable. 

In the case of d = 2 dimensions, we can choose the following orthonormal diagonalizing basis: 



asym 



1 

V2 



1 
-1 



symO 



_ 1 

72 I -1 



^ 1 

u 



seal 



^ 1 

lO 1 



(11) 



In the case of d = 3 dimensions, we can choose the following orthonormal diagonalizing basis: 



asym 



symO 



seal 






2.6 The 2D case with broken symmetries 

It is interesting to study more in detail the case of two-dimensional flows where the velocity field's 
statistics' translation invariance or parity invariance or both are not assumed. We will however still 
assume isotropy. 

As mentioned in Sect. 12 .H^ translation invariance is broken if we don't have the symmetry of Cikji 
under exchange of k with I, or equivalently of i with j. This allows us to have different coefficients 
for 6ikdji and 6ii5jk in Cikji- 

Not supposing parity invariance allows for introduction of an additional term in Cikji, which can 
be written as eikSji + ejidik (recall that simultaneous exchange of i with j and of k with I has to leave 
Cikji invariant). Other expressions of the kind (one e and one 6) give either a zero tensor or the 
same one (up to a possible sign change). In particular (see next paragraph for idea of simple proof) 

^ikSji + ejidik = eii6jk + ejkSu 
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so that this term does not lead to breaking spatial homogeneity of the velocity field. Expressions 
with two e can be expressed with 6s only. 

The covariance matrix Cikji is completely characterized by its action as defined in ifTIH) . and this 
action is linear. In particular it can be expressed as a matrix in the basis of 2 x 2 matrices of (fTTjl . 
For convenience, we order the basis elements as follows: the first two are those in symO , the third is 
of asym and the fourth is that of seal. Then we have the following actions: 



Si j Ski 



/I 0\ 

10 

10 

Vo 1/ 



6ikS. 



ikOjl 



/O 0\ 





Vo 2/ 



6u6 



ilOjk 



(I 
1 


Vo 



0\ 



-1 

1/ 



/O 0\ 



2 

Vo 2 0/ 



This representation is convenient also for finding identities between different expressions. For example 
we have the action 











o\ 




















2 





vo 








0/ 



(12) 



from which we see that 

^ik^ji = SijSki — Sii6jk 

Our next task is to find out when is Cikji, constructed from the above four tensors, a positive 
matrix. For this, we see first that in general (for generic coefficients a, b,c,d£ M) we have the action 



aSijSki + bSikdji + c6u6jk + d{eik8ji + ejiSik) : 



(a + c 

a + c 

V 



\ 



a — c 2d 

2d a + 2h + cj 



The necessary and sufficient condition for Cikji to be positive is that all eigenvalues of the above 
matrix are positive. These eigenvalues are readily found to be 



a + c (twice) , a + h+ {h + cY + Ad^ 

The term under the square-root is always positive, and since the square-root itself is defined to be 
positive, it is enough to retain the two criteria 



a + c > 



a + 6 - \/(6 + c)2 + 4^2 > 
as necessary and sufficient for CjfcjZ to be positive. 



(13) 
(14) 



3 Inertial particles in the linearized flow 
3.1 Basic equation 

Let us turn to the inertial particles. An inertial particle is one whose velocity does not necessarily 
coincide with the speed of the fiuid fiow surrounding the particle. Instead its velocity relaxes to the 
latter by viscous friction, at an exponential rate whose time coefficient is the Stokes time, which we 
shall denote r. 

It is not completely clear what equation to use to describe inertial particles. The work of Maxey 
and Riley is often cited and ad hoc simplifications of their formula are made. Also the fact that 
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describing particle separation is a different matter from describing single-particle motion is mostly 
neglected. Here we shall admit that for very heavy particles (the case of aerosols) the following 
equation is a good approximation both for particle movement and for extracting from it the equation 
on particle separation. Formally, we thus take a system where inertial particle evolution is described 
by the very simple equation: 

§-v. §--',iv-m,^» (15, 

Remark that, in the case of the Kraichnan model, (fTK|l is in fact a stochastic differential equation 
(an SDE; see Sect. l3.4l for a justification of why lfT5|l is an SDE: the idea is illustrated on (|T6|l ). so in 
principle we should also specify what interpretation (Ito , Stratonovich, etc.) we take, but arguments 
of Sect. 13.51 applv here too to show that the particular interpretation convention doesn't matter. 

Eq. lfTK|l gives for the infinitesimal separation r and relative velocity v of two (infinitesimally 
separated) particles: 

dr ^ d-u 1 ,^ , , ^, 

-=., _ = (t).) (16) 

where a{t) is as defined in and in particular depends on the trajectory R{t) of the reference 
particle, so that in principle l[TH|l is not autonomous and should be coupled to (fTK]l. But for the 
purpose of calculating statistical averages of the particle separation it is quite enough to know just 
the distribution of the process cr{t). As was already discussed in Sect. 12. 2( in the case of a Kraichnan 
velocity field the process a{t) is centered Gaussian white noise with two-point function given by (EJ, 
and this is all that we need to know. Once again, for the Kraichnan velocity field, ifTHjl is an SDE, 
and again interpretation convention doesn't matter, as developed in Sect. IH.4I and 13.51 

3.2 Anderson localization form 

The two first-order differential equations of (fT?)]l can be combined to one second-order on r 

d^r 1 dr a ^ 
dt^ " ~Tdt ^ t"" 

or, introducing ipit) = e*''^'^f to eliminate the first-order term: 

d>^ /I a\ y 

Introduce 

E = —1/At'^ (energy), V = a/T (matrix white-noise potential), (17) 

consider t to be position rather than time, and we have the Schrodinger eigenvalue equation associated 
to the Anderson localization problem (though the wave function ^|) takes values in not simply in 
M, and the potential V is accordingly matrix valued, not simply real): 

-^ + V^=EiJ (18) 

This analogy between the system (fT6|l and the Anderson model is hinted at (and used) in [2]. 

In the one-dimensional case (our d = 1), the Lyapunov exponent associated to the growth-rate of 



(more exactly of 

VW+^) is known (cf. Sect.EIHI), but generally not in the multi-dimensional 



case. 
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3.3 Analytical form 

The equation of evolution may also be written for matrices TZ and V instead of particle separations 
r and v: 

dt dt ^ ' ' 

where every corresponding column of TZ and V correspond to a couple f and evolving in the same 
velocity field but independently of the other columns. 

It is straightforward to see that the matrix Z = VTZ^^ verifies the first-order equation 

dZ ~n 1 ~ a 
= -Z^ --Z + - 

dt T T 

since for any matrix-valued function dTZ^^ = —Tl~'^{dTl)Tl^^. We may also introduce Z = Z + l/2r 
and write 

dZ ^2 ^ ^ 

3.4 Passing to an SDE 

Since a{t) is a finite-dimensional Gaussian process, it may be represented as linear combinations 
of a finite number (at most d?) of independent white noises. Denoting the latter by dw (with 
{dwi,dwj) = 5ijdt) clearly we may write, for some (not necessarily square) matrix B{f) 

T-'^a{t)r = B(r)dwldt 

or equivalently 

T-^dSr = B{r)dw (19) 
With this in mind, we may write the block-matrix form equation 



The merit of this formulation is that we immediately see that we have a stochastic differential 
equation. 



3.5 Whatever convention (Ito etc.) 

Furthermore it doesn't matter what convention (Ito , etc.) we use to interpret {SOI, as we shall now 
show. Indeed, write the generic equation 

dX = Adt + Mdw (21) 

where A and B depend on X. Then from passing from Ito to Stratonovich convention, we have the 
additional term \{{{{Mdw) .V)M)dw) , easier to develop in coordinate notation: 

111 

-{Mijdwj{diMki)dwi) = -{MijdiMki){dwjdwi) = -MijdiMkjdt 

since (dwjdwi) = 5jidt. It is easy to check that in our case the above is since Mijdi will only involve 
derivatives with respect to v whereas M^j depends only on r. 

The same reasoning as above applies to single-particle motion described by (|T5|l . It is then 
interesting to note that the r ^ limit may be taken and we thus get a unique convention both 
for the passive tracer particle and the separation of such particles, which can be seen to be the Ito 
convention (cf. Appendix EJ- In other words, in a Kraichnan flow a passive tracer that is the r — > 
limit of an inertial particle is solution of the advection SDE with Ito convention. 
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3.6 The Markov process generated by the SDE 

Now let us look at the Markov process generated by an SDE of form JS^)- We may write the 
probability density function of X as P{x) = {6{X — x)), and use the Ito formula 



df{x) = dx,{dif){x) + - {dx.dx j){didj)f{x) 



to obtain 



dP{x) = d{6{X - x)) = {d6{X - x)) 

= -Aid.,M^ - ^)) + ^{M,kdwkMjidwi)d^^d.,^{S{X - x)) 



-Ai{diP){x) + -MikMjk{didjP){x) 



dt 



since {dwkdwi) = 6kidt. Thus for the Markov process generated by the SDE ipijl it is not M and w 
itself that are important, but only MikMjk, in other words MM"^ , so that lf2Hl may be replaced by 
the effective equation 

dX = Adt + Mdw (22) 

for some driving noise dw, not necessarily of the same dimension as dw, with (dwidwj) = 5ijdt. Note 
how this fixes only the number of rows of M (to be the same as the dimension of X), but not the 
number of its columns. 



3.7 Reducing the dimension of the driving noise 

Applying the general considerations of the preceding section to our case, with the particular form 
Jin of ini), i.e. having 

A=( "'1 M ' » 



V \-\v \B{f) 



we have 

'0 

B(r]B^{r 



MM^ 



and using formula (jHH) for B{f) we get 

{BB^)i.j(f) = Bik{r)Bjk{f) = {BikdwkBjidwi)/dt = (d5,fcrfcd5j7n)/dt (23) 

ID 

= r-^C,k,jinn = ^[{d + l- 2p)r^6,j + 2(pd - l)rirj] (24) 

We may choose any other B{r) with real coefficients giving the same result, i.e. such that: 

Bik{r)Bjk{r) = Bik{f)Bjk{r) = T^'^CikjiVkn (25) 

and use it in by posing 



M 





B{r)^ 

Notice that for any orthogonal basis formed of r and d—\ other vectors of length r, whose matrix 
we shall denote R (i.e. the columns of R are the vectors of the basis), we have 



BB'^ = -^[{d + 1 - 2p)RR^ + 2(pd - l)(i?ei)(i?ei^^i 
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since RR^ = diag^(l, . . . , 1) and r = Rei; then note eief = d\a,g^{l, 0, . . . , 0) and get 

ID 

BB^ = — ifidiagrf((2p + l){d -l),d + l- 2p,...,d + l- 2p)R^ 
(note (i + 1 — 2p + 2{pd — 1) = (2p + l)((i — 1)) so that in particular we may take 

B = R diagrf( v//?Z, • • • , \f^) (26) 

where we have introduced 

/3L = 2r-2z)i(2p + l)(d-l) /3jv = 2T-2L)i(d+l-2p) (27) 

The effective SDE we obtain for is driven by a d-dimensional white-noise instead of the 
d^- dimensional white-noise du;, and is written 

3.8 An even stronger reduction 

Due to symmetries of the problem, a closed evolution equation may be written on the magnitudes 
and relative angle of particle separation r and relative speed v. Let us begin by introducing the 
quantities 

X = v.f Y = f.f Z = v.v (29) 

These verify the evolution equations (note that throughout we use ltd convention to interpret the 
SDEs below) 

dX = dv.f + v.df = [ — v.f + v.v\dt -\ — f.{dSf) 

T T 

dY = 2f.df= 2f.vdt 

2 2JD 2 
dZ = 2v.dv + {dvi,dv^) = [ v.v + -^id + 2){d - l)r.f\dt + -v.idSf) 



where in the last line we used 



^dvi.dvi) = {dSikrk,dSiiri) = {dSik,dSii)rkri 

= Cik,u dtrkri = 2Di[{d + 1 - 2p)d + 2{pd - 1)] (f.f) dt 
= 2Di{d + 2){d-l)Ydt 

It is convenient to introduce the new noises d?/i, drj2 defined by dr/i = f.(d5f) and dri2 = f.(dS'f). 
We only need to know their correlations: 

(dT/i,d?7i) = {dSikrirk,dSjirjri) = {dSik,dSji)rirkrjri 

= Cikji dtnvkrjn = 2Di{2p + l){d - 1) (f f)^ dt 

In the same manner 

(dr/2,d7?2) = {dSikVirk,dSjiVjri) = CikjidtVirkVjri 

= 2Di[{d + 1 - 2p)(f f)(f.f) + 2{pd - l)(f.f)^] dt 

and 

(dr?i,dr?2) = {dSikriVkASjiVjri) = dkjidtrirkVjn 
= 2Di{2p +l){d- 1) (f.f)(f f) dt 
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To sum it up 



dX = [--X + Z]dt+-dT]i 

T T 

dY = 2Xdt 

dz=[ — Z + ^{d + 2){d - l)Y]dt + -dm 



with 



{dm, dm) = 2Di{2p + l){d-l)Y^dt 

(dmAm) = 2I?i [{d + l- 2p)YZ + 2{pd - dt 

(dr/i,d?72) = 2i?i(2p + l)(d-l)Xydt 

We can now go on by introducing 

A = X/Y B = Z/Y (30) 



Then 

dA = d^ = ^dX - ^dY = {--A - 2A^ + B)dt + -d7i 

Y Y Y T T 

y '\ y 9 9 r) 9 

dB = d- = -dZ - -p.dY =\—B- 2AB + -^{d + 2){d - \)\dt + -d72 

Y Y Y ^ T T 

with d7i = Y^^dm and d72 = y~"'^d72, so that 

(d7i,d7i) = 2I)i(2p+l)(d-l)dt 

(d72, d72) = 2Di [{d+l- 2p)B + 2{pd - l)A^] dt 

(d7i,d72) = 2Z?i(2p+l)(d-l)Adi 

Introduce C = B - A^ . Then 

dC = dB- 2AdA - (dA, dA) = [ — C- AAC + -^{d - l)(d + 1 - 2p)l dt + -d73 

r T 

where we have introduced d73 = d72 — Ad7i. We have 



(d73,d73) = (d72,d72) + yl^(d7i, d7i) - 2A(d7i,d72) 

= 2Di{d+l - 2p)Cdt 
(d7i,d73) = 



Finally introducing F = C^/^ we arrive at 



dF = dCi/2 ^ ic;-i/2dc - -C^y^{dC,dC) 
2 8 



[-(- + 2A)F + ^{d-2)id + l- 2p)F-^] dt + -d74 



with d74 = C ^/^d73 so that 

(d74, d74) = 2Di{d+l- 2p) dt 

At this point, let us introduce z = A + iF + 1/{2t) (so the real part of z is the (signed) component 
of V along r, divided by r and then shifted by l/(2r), while the imaginary part of z is the norm of 
the component of v perpendicular to r, divided by r) and d7 = d7i + idj^ to have 

dz = [^-z^ + i^{d - 2){d + 1 - 2p)(9z)-i] dt + id7 
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with (since (d7i,d74) = 0) 



(d7,d7) = (d7i,d7i) - (d74,d74) = 4L»i(p(i- 1) dt 
Or, using definitions from lfT7|l and l(?7|) . we can write 

dz= [-E-z2 + i^^^/3^($5z)-i]d* + d7 (31) 

with 

d7 = d7L + i djN (32) 
where dji and d7Ar are real white noises with correlation 

(d7L, d7L) = (3l dt {d^N,dlN) = dt {d^iAlN) = (33) 

so that in particular (d7, d7) = (/3l — dt. We emphasize that 7 and 7 are Brownian motions 
not on the real line but on the complex plane, each equivalent to a two-dimensional real Brownian 
motion. Formula l(3T1) (with other notations) appears in jH]. 

We then see that the case d = 2 is special since then the equation on z is analytical. For d > 2 
the additional drift comes from the fact that we reduced the diffusion of v perpendicularly to r to its 
radial part, exploiting rotational symmetry around the direction of r. Then we have the well known 
inverse-radius drift term of the radial Laplacian. 

Also the case {d + 1 — 2p) = is special, but as we will see in the next subsection, it is a much 
more particular case, since v tends to align with f so that we recover the one-dimensional situation 
which is explicitly solvable. 

Since we also have (cf. Eqs. pUl l 

dloglrl = — dY = Adt = (^z - —) dt 
^' ' 2y V 2r/ 

we have the expression for the Lyapunov exponent of particle separations 

-(^)-<«^>-^ ,34, 



3.9 A solvable case 

The "physical" values for [3n / are those that correspond to < p < 1, i.e. | < (3^/(31 < 1 + ^zjj 
in the sense that only these correspond to the evolution of the separation of inertial particles in a 
linearized Kraichnan flow. However, the reduced SDE (|28|l or pTjl is meaningful as long as /Sat > 
(cf. respectively ^ and ((SSI))- 

The case for which the Lyapunov exponent is known is when jS^ = 0, and additionally the initial 
condition is such that v \\ r, as then (pTTjl has real initial z and both the imaginary term in the drift 
and in the diffusion (cf. ((32|l . (jSSl) disappear. We thus recover the classical Anderson localization 
problem: introduce il) = e*/^'^r (here r E M is the signed coordinate of r along the line defined by the 
initial condition r(0)) verifying 

, d7 

-d^+d^^ = ^^ 

then z = ip' /ip verifies H3H) . so {z) is the growth rate of InlV'l which we know from the Anderson 
problem. Hence, recalling the simple relationship ((M|l between (z) and A, the latter may be expressed 
in terms of the Airy functions of first and second kind, Ai and Bi respectively, and their derivatives 
Ai' and Bi' , as 

1 



' 2r 



_^^^ Ai'{c)Ai{c)+Bi'ic)Bi{c) 
Ai^{c)+Bi^{c) 



(35) 



where c= l/[4T2(/3i/2)2/3]. 
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Later on we need to examine if the v \\ r initial condition is really necessary for the above formula 
to hold. Numerical simulations of Ijlftp (in fact of (jl]8|l . which is equivalent even in d ^ 2 dimensions, 
since (fT6|) reduces to pTj) . which for = is identical to and the latter is derived from l(38|l 
and is sufficient to give the corresponding Lyapunov exponent) show that (at least for the examined 
values of parameters), this condition is irrelevant. 

It is tempting then to think that this irrelevance of the initial alignment of v with f is due to 
the asymptotic "collapse" of v onto the direction of r at long times. In other words, that z, even 
if initially it has an imaginary part, can be considered to be real at large times, because it would 
be "attracted" to the real axis by the dynamics of its evolution. However this is not the 
evidenced by numerical simulations, since \?R.z\ has a finite limit when t — > oo, which in fact we can 
also calculate. The more general treatment of Sect.Elwill shed light on this and justify the expression 
l(35)l even for non-aligned initial v and r. 



4 The 2D case 

4.1 Passing to the complex notation in 2D 

Recalling Sect. IH.7l and in particular the requirement (see also the last equality of l(2U), we can 
choose for B{f), in the 2 dimensional case, the particular form 



PLr2 



with 



Pn = 2t~^Di{3-2p) 



(36) 

This case is special, since R (the matrix of some orthogonal basis whose first vector is r; cf. Sect. 13.7(1 
can be given as a linear function of r, namely we chose above R = (r, 0^/2^); where stands 
for rotation by angle 7r/2 in the positive direction. In general there is no such linear representation 
(however for dimensions d that are a power of 2, the Cayley-Hamilton algebra of dimension d permits 
to do a similar rewriting). Then 



B(r)dw 
Thus, if we define 

so that 



'T3^dw2 ^/JJJ^dwl J \r2 



a 



fJJldwi -y/PNdW2\ 

^dw2 y/PEdwi J 



afdt = B(r) dw 



then can be written as the system 

dr 

dt ~~ 



dv 
dt 



1 ^ 

■—V + ar 

T 



(37) 



(38) 



On the other hand we see that, if we identify to C in the standard way, then a corresponds to 
multiplication by (and we re-employ a by abuse of notation) 



(39) 



We may pass altogether to the complex notation 



dr du 1 _ 
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where r and v are the complex equivalents of r and v respectively. 

We may proceed further analogously to Sect. 18 .21 keeping E from lfT7)l and defining V = a/r 
which is now a complex-valued random potential. Then ili = e*/^'^r verifies 



+ V'il) = Eip 



(40) 



similar to (|TH|l . 

On the other hand one may proceed along the lines of Sect. \',\.\\\ introducing z = v/r. Then z 



verifies dz/dt 



z/t + a and 2 = 5 + l/2r verifies 
dz 9 1 



(41) 



Note finally (notwithstanding that r is complex) 

A = (Aiogr) = (^) 
dt r 



2r 



+ (-) 



4.2 An other solvable case 

Following ideas of jSj, we will now show that in the 2D case the Lyapunov exponent may be explicitly 
computed also when /3i = 0. The important particularity of the Pl = (i-e. p = —1/2) situation 
is that the steady-state distribution of z is concentrated entirely on the complex half-plane with 
diz > 1/(2t), due to the fact that on the line ?R.z = l/(2r) the velocity field always points inside the 
above-mentioned half-plane. Indeed, writing z = x + iy = 1/(2t) -|- iy we have the velocity field 



-E 



1 



{2Tf 



2r' 



and > except for y = 0. 

In this particular situation the Laplace transform F{p) = (exp(— pz)) is well defined for p G M+. 
Using the Ito formula 



de" 



-pz 



-pe-P^dz + ^e~P'{dz,dz) 



and substituting dz according to (02), we get for the steady state 



dt 



Fip) 



{E + z^)p- 



Pn 2 



-pz 



P 



dl + E 



-p 



F{p) 



Since the above is defined only for p > 0, we may safely divide by p and retain simply 







dp+E- —p 



F{p) 



At p = we have the obvious boundary condition -F(O) = 1. For p going to -|-oo we need 
limsup|F(p)| < 1 since in the steady state |exp(— pz)| < 1 almost surely for p > (due to z being 
almost surely in the positive half-plane). These two boundary conditions are in general sufficient to 
uniquely determine F{p) in the steady state. 

This differential equation can be readily solved (e.g. by introducing u = p — 2E/ (3^ it reduces to 
the Airy equation (5^ — {j3N /2)u)f = 0). The solution that decays (in fact the only one that does 
not grow exponentially) when p — > -|-oo and that verifies F(0) = 1 is 



F{p) 



Ai{-{ 



2 . 



Ai{-m-V^E) 
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where Ai stands for the Airy function of first kind. Finally we arrive at 



{z) = -{drm 



Introducing c = l/[4T^(/?Ar/2)^/^] = -(^) '^^^E (note that preceding definition of c was with 
instead of Pn^-), we have for the Lyapunov exponent 



^_ ^-1/2 Ai'(c) 



Ai{c) 



(42) 



4.3 Decay of Lyapunov exponent at large Stokes numbers 

Let us introduce the adimensionalized Stokes time 

St = DiT 

which is nothing else than the Stokes number for the case of a smooth Kraichnan velocity field. We 
are interested here by the decay of the adimensionalized Lyapunov exponent X/Di when the Stokes 
number St goes to infinity. 

Numerical results (see subsection 16 . II and figure ISJ seem to show us that the Lyapunov exponent 
A of pair separation varies monotonically with the compressibility degree p of the advecting flow, 
notably the larger p (i.e. the more compressible the flow) the smaller is A. Thus the two extreme 
cases of p = —1/2 and p = 3/2 bound the other, intermediate cases. 

It is then noteworthy that for large r the exponents given by both (jH^ and decay at the 
same rate. Indeed, we have Pl,Pn ~ Dit~^ so in both cases c ~ St'"^/^, and since Ai{c), Bi{c), 
Ai'ic) and Ai' {c)Ai{c) + Bi' {c)Bi{c) all have finite non-zero limits when c — > 0, we have the behavior 
\/Di ~ St-'^/^ in both cases, so this should be also the general situation. 

An other possible way to arrive at this result is by taking r to oo while fixing a in dHJ. At the 
limit, when 1/r = 0, we have simply z = z and ll^Tj) reduces to 

dz 9 
dt=-^ 

It would then be enough to see that this evolution equation leads to a (unique) stationary state with 
{z) finite and non-zero, and that this is indeed the r — > oo limit of (z) (i.e. commutation of t — > oo 
and r — > oo limits). 

Indeed, since a is fixed, this means that in fact I3l,I3n are fixed, so that (recall JSEl) Di ~ 
and St ^ T^. Since we suppose that for fixed a the Lyapunov exponent A goes to a constant, the 
adimensionalized Lyapunov exponent X/Di behaves like D^^ ~ ~ St~'^^^: 

— ~ st-'/^ 

Di St^oo 

5 Other solvable 2D flows 

The idea in this section is that the above two special cases of solvable situations can be slightly 
generalized. The cases we could solve were when a of was either y/j3Ldwi/dt or iy/ (i^ dw2/dt. 
In fact we will now deal with the more general case of 



a = V/3dw}/dt (43) 

for /3 € C (and dw is the standard white noise). 

By analogy with Sect.^ where + = 8-Di/r^ (from ((SEl)), we define here the characteristic 
inverse time-scale of the flow as 

Di = -1/31 
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5.1 Solvable flows 

Let us start by finding the 2D flows which lead to such a. Noticing that one has (analogously to 

Bik{f)Bjk{f) = {BikdwkBjid'Wi)/dt 
and using to re-express the left hand side and l(H7jl to re-express the right hand side, we get 

T~'^Cikjirkri = {aikrkdt, djindt) /dt (44) 

so if we deflne 

Cikji = {d'ikdt,d-jidt) /dt 

then 1^ can be written 

T~'^Cikjirkri = Cikjirm (45) 

Considering the complex a of l|i3|l as a real linear transformation acting on C = M^, its matrix 
is written (also a by abuse of notation), with help of the notation ^/P = a = ai + iajy [a £ C, 
ttL, UN G K) 

a=f«^ -""Adw/dt 

In coordinate notation this gives 

^ik = {oiL^ik - ttNeik) dw/dt 

Leading to 

Cikji = {oLL^ik - aNeik){aL5ji - UNeji) 
and using (Tra . this further gives 

= al6ik6ji + a%{6ijSki - SuSjk) - aiaNieikSji + ejiSik) 

With the above form for Cikji and using the most general form respecting isotropy for Cikji (cf. 
Sect. ESI), that is: 

Cik,ji = a5ij5ki + hSikSji + cSuSjk + d{eikSji + eji6ik) 
with arbitrary a, b,c,d£ M, equation becomes 

T~'^{a6ijr'^ + (6 + c)rjrj + d{eikrjrk + ejkrirk)) = 

a%Sijr'^ + («L - ot^^i'^j - OLLaNieikTjrk + ejkriVk) 

Coefficients are identified 

a = 6 + c = a| — a% d = —aic^N 

and since only the sum of b and c is determined, let us write b = aj^ — s, c = s — aj^ for some s G M. 
The positivity conditions ifT^ . lfn|l give 

s = a + c> 
-s = a + b- 7(6 + c)2 + 4^2 > 

meaning that the only possibility is to take s = giving 

Cik,ji = [alSikSji + a%{6ij6ki - Sudjk) - aLaN{eik^ji + eji^ik)] 
= T^iaiSik - aNeik){aLSji - aNeji) 

Such Cik,ji can occur only in a flow (cf. Sect. 12.61) where parity invariance is broken because eikSji + 
Sji^ik appears, and spatial homogeneity is broken because 6ikSji and 6ii6jk have different coefficients 
(namely r^a^ and — r^a^, which are equal only in the case tol = ton = 0, i.e. Cikji = 0: a very 
trivial case). 
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5.2 Calculation of the Lyapunov exponent 

This is very much hke in Sect. 14.21 




Figure 1: The complex plane of z is represented. There is diffusion only in a single direction, parallelly 
to a = \f^. The dynamics of 2;, governed by (|lT|l with a as in (|^ . is such that z can cross the boundary 
line betvi^een the hashed half-plane (unstable) and the clear half-plane (stable) only moving out from the 
hashed half-plane and into the clear one, whereby asymptotically the distribution of z will be supported 
entirely on the stable half-plane. EI_ is the lower half-plane. 

Let us first establish the existence of a stable half-plane for z, by which we mean one on whose 
boundary line the drift velocity of z points always towards the interior of the half-plane. We may 
suppose, without loss of generality, that < arg/J < vr and in consequence < arga < 7r/2. We 
then claim that the half-plane 

z-^ 

{z : ^ G H„} 

a 

with ]HI_ denoting the complex half-plane with negative imaginary part is such. Note that this 
half-plane may be parametrised as 

— + fia - iua fi eR+ (46) 

2r 

In particular its boundary line is parametrised as l/2r -|- fj,a with /i G M, and —ia is perpendicular 
to the the boundary line and points to the interior of the stable half-plane. 

What we need to show is that, on the boundary line, the scalar product of the drift of z with —ia 
(considered as a vector of M^) is positive. Remark that, viewed as vectors of M^, the scalar product 
of two complex numbers x and y may be written as ^{xy). Then we have for z = l/2r -|- fia 



-E - z^){-ia)] = 5R[( ^la - ^i'a'){-ia)] 

T 



/iz|a|^ — ^^|a|^ia] = //^|apK(— za) 

r 



The idea is again that in the stationary state the distribution of z will be supported by the stable 
half-plane, being identically zero outside of it. In this case the Laplace transform F{p) = (exp(pz)) 
is well defined for p € (ia)~^M+ since then K(pz) < 5?(l/2rp) as can be easily read off from 
making | exp(pz)| uniformly bounded in z for fixed p. 

As previously, in the steady state we may write 



o = Pip) 



-{E + z^)p+^p^ 



: p{-dl -E + ^p)F{p) 



We may divide the above by p without loss of information, since at p = we have the boundary 
condition -F(O) = 1, and as \p\ — > 00 with p G (^0)"-*^]^+, the F(j)) has to decay, and these two 
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conditions determine uniquely the solution. 

id'^+E-^p)F{p) = (47) 
This equation is related to the Airy equation, if we introduce the new variable u = p — 2E/(3: 

{dl - ^u)f{u) = 

of which two independent solutions are (formula 10.4.1 in [T|) Ai{{l3 /2Y/'^u) and 74i(e^'^*/'^(/5/2)^/^ti), 
where Ai denotes the Airy function of first kind, and {(3/2)^/^ is taken with the standard definition 
of the root, (i.e. real positive on the positive real axis and having branch cut along the negative real 
axis). 

For large argument, to lowest order, the asymptotic expansion of the Airy function Ai reads 
(formulae 10.4.59 of [T]) 

^^(0 = + 0(C-'/')) (I arg CI < vr) (48) 

Let us apply the above with C = {l3/2fl^u = {f3/2f/^{p-2E/f3). SincepG (ia)-^M+ and a = ^/^ 
(with standard definition of the root), when |p| — > oo we have argn — > argp = — (arg/3)/2 — 7r/2, 
whence argC — > — (arg/?)/6 — 7r/2, so that for large enough \p\ we have — 7r/2 > arg^ > — 27r/3 (in 
particular | argC| < vr), and — 37r/4 > arg(C^/^) > — vr. This implies, when \p\ oo, that 3f?(C^/^) < 0, 
and since |C| ~ {\f3\/2)^/^\p\, we get from ^ that \Ai{{P/2y/^{p - 2E/(5))\ oo. 

Repeating the same treatment with C = e^''^/^{(3/2y/''^u = e^'^'/^{(3/2)^/^{p - 2E/ (3), when \p\ 
oo we have argC — > — (arg/3)/6 + vr/G, so that for large enough |p| we have < arg^ < 7r/6 (in 
particular |argC| < vr), and < arg(C'^/^) < 7r/4. This implies, when |p| oo, that 3f?(C^/2) y 0, 
and since |C| ~ (|/3|/2)V3|p| ^ ^^^^ ggj ^j^^^^ \Ai{e^'''/^{(3/2)^/^{p - 2E/P))\ 0. 

Hence the unique solution of l|17|) verifying the appropriate boundary conditions at and at 
infinity is 

Aiie^-^/H(3/2)y'{p-2E/(3)) 
Af(e2W3(/3/2)i/3(_2£;//?)) 

Remark that {z) = 5pF(p)|p=o, giving 

iz) - e2.V3.^/,.i/3 ^^-(e^"/^(/3/2)^/^(-2i^//3)) 

^ ' Ai(e2W3(/3/2)V3(_2i?//3)) ^^^^ 

We see (if we evaluate the above formula for non-real /3), that {z) is in general complex valued. By 
analogy with Sect. 0] we define 

The real part of A once again gives the exponential growth rate of the pair separation r, while the 
imaginary part of A gives the average rotation speed of r. In particular this mean angular velocity is 
in general non-zero. 

Eq. (gni) may be further simplified by noticing that e^'^'/^ifi /2Y)^ {-2E / (3) = -Ee^'^'/^{I3 /2)-'^/'^ , 
and |e2W3(^/2)-2/3| = |/5/2|-2/3 and arg(e2W3(^/2)-2/3) = 27r/3 - (2/3) arg /? = (2/3)(7r - arg/3), 
and since < arg /5 < vr we have vr— arg/3 = arg(— 1//3) = aig{{—j3)^^) so that arg(e2'^*/^(/3/2)~2/3^ = 
arg((-/3)-2/3)_ From this we conclude that e2'^*/3(/3/2)-2/3 = (_/3/2)-2/3_ Along the same lines, 
arg(e2'^*/'^(/3/2)^/^) = 27r/3 + (arg/3)/3 = (arg/3 — 7r)/3 + vr. However we have to pay attention here 
since arg/3 — vr = arg(— /3) is true only for < arg/3 < vr but not for arg/3 = 0, so to sum it up 
arg/3 — vr = arg(— /3) is true only for (5 ^ M+. Thus, except for the case of (5 real positive we have 
(arg/3 - vr)/3 + vr = arg(-(-/3)V3) go that e2W3(^/2)i/3 = _(_/3/2)i/3. Thus we also have 



19 



Note that throughout we have supposed f3 to be in the upper complex half-plane. However this last 
formula has the interesting feature that it is real, so that substituting /? with its complex conjugate 
/3, we also get for (z) the conjugate value, and from simple symmetry considerations with respect to 
the real axis, this is the correct answer. Thus formula (|5?Hl applies for all (3 £ C \ ]R+. 

To understand the problem with (3 € M+, note that in this particular case, when the diffusion is 
parallel to the real axis, there are two stable half-planes: the lower and the upper complex half-planes. 
The equilibrium distribution of z depends on the initial distribution of z, since z cannot cross the 
real line. As we have already showed, for z starting from the lower half-plane (this was only implicit 
in our considerations, through the fact that we considered the stable half-plane to be the lower one) 
we arrive at formula l|49p . Starting from the upper half-plane we arrive at its complex conjugate. 
Hence the natural branch cut in l(5?Hl for f3 S M_|_. 

We recover the two solvable cases presented in Sect. 13.91 and Sect. 14.21 The former one, with 
Pn = (hence diffusion only in the real direction), corresponds to /? = /3l £ R+. Using (|l9|l and the 
identity (formula 10.4.9 of [T]) 

leading also to the identity 

^•/(-g27rV3^) = e-'^^^I^^Ai(e^''^I^C) = -e-^'/'\Ai'{C) - iBi'(C)) 



we have 



Ai{-E{pL/2)-^'^) - iBi{-E{l3L/2)-V-^) 

Introduce c = —E{(3l/2)~'^/^ = 1/[4t^(/?l/2)^/^], and multiply numerator and denominator in the 
above displayed formula by the complex conjugate of the denominator, finally use the Wronskian 
identity (formula 10.4.10 of see Index of Notations p. 1045 of the cited source for sign convention 
for the Wronskian) 

Ai{0Bi'{0-A^{QBi{(:)=7,-^ 



in order to get 



(z) = —c 



1 ^_^/^Ai'{c)Ai{c) + Bi'{c)Bi{c) - iir"^ 



2t Ai^{c) + Bi^{c) 

So what is the meaning of this imaginary part ? Recall that our initial assumption was that z 
is inside the stable half-plane at asymptotically long times, and as a limiting case of arga \ this 
half-plane is that of negative imaginary parts. The asymptotic distribution of z is indeed supported in 
this half-plane if the initial z is in it. If however we start with z initially distributed with probability 
1/2 on each side of the real axis, or if we add infinitesimally small diffusion in the imaginary direction, 
then we get in the upper and lower half-planes distributions which are mirror images of one-another, 
and in the lower half-plane it is 1/2 times that of the case when there was nothing in the upper 
half-plane. 

That (z) does have the imaginary part predicted above can be readily verified in simulations (see 
Fig. El the f3 = 1 curve), even if we use the scheme of infinitesimal diffusion in the direction of the 
imaginary axis, by evaluating and confirming the relation 

{\QZ\) = ic-V2. 



2t At^{c) + Bi^{c) 

We also recover the case /3l = (hence diffusion only in the imaginary direction) corresponding 
to P = —(3n G Here we use (|Kn|) and get 

(A- (R /ou/3Mz^(M^):!^!) 
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or, introducing c = -E(/3jv/2)~2/^ = 1/[At'^{I3n 



^ > ~ 2r Ai{c) 



6 Numerical results 



In this section we focus mainly on the d = 2 dimensional case. 



6.1 Naive simulation of inertial particle separation 

In order to get the Lyapunov exponent of pair separation, the most straightforward thing to do 
is to simulate the reduced equation l(88|) . Note that this is better than simulating since at 

moments when particles nearly collide (happens when particles are moving towards each other with 
small impact parameter (compared to ru)), the quotient v/r = z — 1/2t turns very fast, which would 
require additional handling in the code, whereas in fact v evolves little at near collision since v r jr 
and T evolves linearly for which no special care is needed in the simulations. 

The convergence of the Lyapunov exponent is quite slow, one can hope for a 0.1% relative precision 
at most, which is not enough to try to guess a functional form through Plouff's inverter, but more 
than enough to disprove some claims of formulae. 

Note also that for r small the system is stiff and convergence gets worse. The equation of evolution 
should be solved analytically between two time-steps, with taking for r on the right-hand side its value 
at the beginning of the time interval. Indeed in this case we have a linear (multi-dimensional) equation 
with time-independent homogeneous term (only the inhomogeneous term contains the driving noise) 
which can be explicitly solved, and we get an effective driving term whose statistics can be computed. 
This is treated in subsection 16.41 

Numerical simulations were done in series with fixed /3l, (i^ (with additionally fixing Pl+Pn = 1; 
note the identity Pl + Pn = 8Di/t'^ leading to = with our choice) instead of fixed Di (note 
that p is fixed by the quotient of Pn)- However we are only interested in the the adimensionalized 
Lyapunov exponent X/Di in function of the Stokes number St = Dit. Simulations were done for 
values of r between 0.1 and 4, which, according to the above, correspond to Stokes numbers between 
1.25-10"^ and 8. Also, for reasons of better numerical stability, instead of taking Pn = 0, we used in 
that case Pn = 10~^. Total time of simulation was 5-10^ with time-step 10~^. 

The plot of numerical results, at the end of the paper, has been cut up into three pieces ([0,0.1], 
[0.1, 1] and [1, 7.5] that are shown on different scales). The vertical axis is X/Di the adimensionalized 
Lyapunov exponent. The curves are ordered in the same top to bottom order as the keys. Solid line 
curves are the analytical solutions, dashed lines are just linear interpolation of numerical data. Data 
itself is represented by boxes, either open or closed. 



6.2 Estimates of v and of the change of r during a time interval At 

Let us estimate first v and then the jump in r during some time interval At. We suppose |r| « 1 at 
the beginning of the time interval, without loss of generality, since the system under study is scale 
invariant. The case we are interested in is when At is the time-step of our simulation, but here 
we talk about the "real" v and r, not the simulated ones. Hence we only consider the case where 
At<^ 

The forcing of is a white noise with intensity proportional to r^^-y/I^ and the "memory" of v 
is of order r, so that ~ T~^^/Diy/T = ^/Di/t. 

As for r, we need to distinguish the cases when r < At and when r > At. The time derivative 
of r is u ~ Di/t. Now if r < At then we can say that v is correlated over a time r so that during 
each interval of time of length r, the vector r changes by approximately ry^ Di/t = \JD\t. There 
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are At/r such independent changes which gives a total change of ^D\T\J IS.t/T = \JDiIS.t, just as if 
r underwent a stochastic diffusion with diffusion constant not depending on r. 

In the other case of r > At the velocity is correlated over times larger than At so that the change 
of r is approximately \v\^t ~ sj Di/rAt. For later use, we remark that this latter can be estimated, 
due to the hypothesis r > At, as y^Di/rAt < y^Di/ AtAt = y/DiAt, so that in this case also r 
changes by at most ^/DiAt. 

6.3 Integrating dlogr^ 

To obtain numerically the Lyapunov exponent A, we essentially make use of the formula 

/ d log r2 \ 

However, as stated above, the change Ar of f during some time interval At can be bounded, in 
general, only as Ar < y/DiAt (but at least this bound is uniform in r). This is fine, but it means 
that to calculate the change of a non-linear function / of r we have to develop / to second order in 
r, which means in other words that we need to use the Ito formula. 

In practice, we have 

dlogr^ = ^2ri dri+ - ^riVj + ^^ij) drj drj (51) 

and this is the formula to use in simulations to get a correct result (in the simulations it appeared 
clearly that for r < At the first-order formula gives completely false results). 



6.4 Overcoming stiffness 

For small values of the Stokes time r the system has a large parameter (v.g. r^^), so that it is stiff. 
In principle one needs to take a time-step At to integrate it that is smaller by one or several orders 
of magnitude than r. However in the r — > limit the resulting process for r is still well defined, as a 
diffusion process, which may still be integrated with a finite time-step At which only has to be small 
compared to the inverse diffusion coefficient D^^ . Can we devise an integration scheme that doesn't 
need At — > as r — > and which is correct both for r < At and for r > At, given At ^ D^^ ? 

The answer is affirmative as we now expose. Since in the r ^ limit the "fast" variable is v we 
will integrate its evolution over one time-step say between and t = At (for ease of notation we will 
use here t instead of the more clumsy At), analytically. This is possible if we keep in the evolution 
equation f fixed to its value r(0) at the beginning of the time-step, since then the system becomes 
linear with a constant evolution matrix and some time-dependent forcing: 




See Appendix for justification of this approximation. 

Let us start out with the study of a generic system of this form 

dX = AXdt + dBt (53) 

Its solution is ^ 

X(t) = e*^X(0) + / e(*^^)^d4 (54) 
Jo 

We thus have a deterministic part and a stochastic one which is linear in the increments of the 
Brownian driving process Bg, in other words the stochastic part is in the first chaos of Bg, hence it is 
a mean zero Gaussian random variable. Also since on (sub-intervals of) different time-step intervals 
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the increments of the driving Bs are independent, we really need to know only the covariance of the 
stochastic part on any given time-step: 



Let us call K{s) = exp(sA), and notice that for a vector n we have the equivalent writings u®u = uif" . 
Introduce the equal-time correlation matrix x such that {dBg dBg') = x5(s — s') dsds'. Then the 
above displayed formula simplifies to 

ft ft 

Kit - s){dBs dB,,)K^{t - s) = (55) 

(56) 



Jo 



[ [ K{t-s)xK^{t-s')6{s-s')dsds' = [ K {t - s)xK^ {t - s) ds 
Jo Jo Jo 



At this point one needs to write the above correlation matrix as UU'^ for some matrix U. This is 
most easily done by way of the Cholesky decomposition which is simply the prescription of U being 
triangular (lower or upper), because the coefficients of U can then be recursively computed by using 
only algebraic operations and square roots. This is algorithmically much simpler than if for example 
we tried to diagonalize the correlation matrix to take its symmetric square-root. 
In our case it is easily computed that 

K{t -s) = 

but note that this a block-matrix notation where each block is proportional to the 2x2 identity 
matrix. As for Xi we can see it also as a 2 x 2 block-matrix of 2 x 2 matrices, whose only non-zero 
block is the lower-right one, and in fact 






^ \0 B{f{0))B^{f{0)) 

Since every block of K{s) (each being a scalar matrix) commutes with the non-zero block of x, 
we can start by evaluating 

L{t,T)= l"K{t-s)(l ^\K^{t-s)ds 



r2[t-i(l-e-7)(3-e-^)] W-e-^)' 



^2 



and its Cholesky decomposition can be given as 



T^/t — 2rtanh^ r(l — e T)^/-|tanh^ 

U ' 

\/i(l-e-^^) 
Finally with the full Cholesky decomposition the driving noise process we need can be generated 

as 

with the r] independent Gaussian white noises, {dijidrjj) = 5ijdt. It is interesting to notice that this is 
tantamount to using to mdependeni realizations of the Kraichnan field {B{dr]i,dr]2) and B{dri3,dr]i)) 
to drive the inertial particle process. It is not difficult to see that this is the general case (i.e. more 
than two particles in a non-linear velocity field may be adequately driven in a finite time-step scheme 
by two independent realizations of the fluid velocity field, just as here). 
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6.5 Extrapolation from moments of positive even order 

Another possibility is to compute the exponential growth rates of positive even order moments of 
the particle separation, which can be obtained analytically (though implicitly only, as the largest 
real root of some polynomial), as described in Appendix^! Then, the curve which we know only for 
even natural numbers has to be extrapolated in order to guess its derivative at 0. This is a tough 
business, and it is hard to get anything better than 10% accuracy. 

There are several ways of calculating the top eigenvalue. For small order 2n of the moments 
(until n PS 10), it is possible to calculate directly the characteristic polynomial of the matrix and 
solve numerically for the largest real root. 

For larger n (up to n ^ 100) it is more convenient to iterate the system x i— > Mx where M is the 
matrix of the system. Since M is a sparse matrix, this is fast and efficient and doesn't require a lot 
of memory. 

For n even larger it is better to use the more memory intensive method of iterating A\-^ with 
Aq = M, since this way the power of M computed is in 2^ for N iterations, whereas for the previous 
method it was just N. 

7 On some conjectured formulae in 2D 

In ^01 appeared the following conjectured formula for the top Lyapunov exponent of inertial particles 
in a 2D Kraichnan flow as described above: 



2^ 



3f? 



1 



2rL»i 



1 + C-V2 



Ai'jc) 
Ai{c) 



(57) 



with 



[4r2(/?^ - /3i)2/3]-i = [16(1 - 2p)5i]-2/3 



(58) 



This would mean that X/Di changes with p only through simple rescaling. Indeed, introducing the 
function 



Fix) = ^ 



1 + X3 



Ai'{x~3) 

2 

Ai{x~3) 



formula (|57|l is equivalent to 



— = (l-2p)F(16(l-2p)5i) 
^1 



(59) 



with the only subtlety that F behaves differently for positive and negative arguments. In fact, we 
have for x > 

,Ai'{x~3) 



F{x) 



1 + X3 ■ 



Ai{x~Ti) 



but 



F(-x) = - 

X 



2 2 2 2 ' 

1 Ai'(x~3)^i(x~3 ) + Bi'(x~3)Bi(x~3) 
1 — X3 ^ — ^ ^ ^ ^ ^ - 

Ai'^{x~3) + Bi'^{x~3) 



To see that (|57|l cannot agree with numerical results, the easiest is to see that (|59|l would imply 
A = for p = 1/2 for all values of St, since F{0) = 2, and clearly this is not the case. 

We compare the conjectured formula with effective values. Notice that for p = the match is 
not too bad, but for larger values of p it becomes worse. 

Later, in jH], it was recognized that (|H7)l cannot be always correct, but it was hypothetized that 
perhaps it could be exact for the p = 1/6 case. To help compare their results with ours, we give the 
expressions of their parameters F, e^, 7 and z in function of our variables: 



3-2p 
l + 2p 



{l + 2p)St 



7 



1 

T 
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where c is that of formula (|58|l . 

They also give a modified formula, which in our notations is 

A = ^(p)e-6(T+W + (1 - 2p)F(16(l - 2p)St) 

where x(p) is an unknown prefactor, that according to the authors of [Hj does not depend on St. 
But they only claim that the additional term should be some sort of leading order correction. 

8 Conclusions 

We succeeded, for certain time decorrelated flows, to find an analytical expression for the Lyapunov 
exponent of inertial particle pair separation. Other cases had to be dealt with numerically. It is 
however interesting that the analytically tractable cases bear lot of resemblance with the general 
case and thus have a predictive power. For example the D'^^X ^ St""^/^ behavior was first noticed 
for the analytically available solution. 

The finding of slightly generalized cases of the Anderson localization problem which permit exact 
solution could be also of broader interest. 

We have described so far only the Lyapunov exponent, but a more general object, the distribution 
(in fact large deviations) of finite-time Lyapunov exponents, contains further useful details of particle 
pair separation. Our work can be extended, following ideas of |TT|, to gain more information on this 
object. 

Hence, the analogy between inertial particles in a random flow and the Anderson localization 
problem (which itself is related to the Edwards polymer model also known as the weakly self-avoiding 
Brownian motion (or random walk)) is fruitful and holds still more promises. 

A Constant r approximation 

Instead of using the approximated equation ifK^ , we can start out from the true system 




which has the generic form (analogously to (EBI)) 

dX = AXdt + dBt{X) 

where it is in fact the covariance oidBt that depends on X{t) as {dBs^dBgi) = x{X{s))5{s — s') dsds' 
What is of interest to us in the numerical simulation is the mean and covariance of 

AX = [ dX = X{t) - X{0) 
Jo 

We can again write a formula of type l(5^ 

X{t) = e'^X{0) + re(*"^)^d4(X(s)) 
Jo 

but now this is an implicit formula since it contains X{s) for s > on the right hand side also. 
For the mean displacement we still have immediately 

{AX) = (e*^ - l)X(O) 
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For the covariance however we only have, instead of (|56|l . the semi-exphcit 

{[AX - {AX)] [AX - {AX)]) = [ K{t- s)x{X{s))K^{t- s)ds 

Jo 

It is at this point that we use the estimates of subsection l6.2( notably that X{s) — X{0) < ^/Dis < 
y/Dii if we take |^(0)| 1. Then we have a relative error by replacing xi^is)) with xi^i^)) of 
order bounded by \JD\t. This basically means that in our approximate method we simulate a 
system whose diffusion coefficient would fluctuate around the real one, with a relative precision of 
order \/D\t. Since in the simulations that are reported here we took D\ = and a time-step of 
t = 10~^, we have y/Diir ~ 3.5 10~^r, and since we used only r < 4, we have \JDitT < 1.5-10"^. 



B The T ^ limit 

When the Stokes time r goes to 0, the movement of an inertial particle solving lfT5|l tends to that of 
a passive tracer that follows the velocity field, solving the SDE 

dR{t) = dUt{R{t)) (60) 

where we preferred writing df/j instead of U to clearly mark that the velocity field is a white noise 
in time. The solution of this equation depends, in general, on the convention we adopt (Ito or 
Stratonovich or other), except if the flow is incompressible or its statistics is locally isotropic. 

However, as pointed out in Sect. 13. 5( for finite r the evolution of the inertial particle does not 
depend on the adopted convention: we don't have any degree of freedom. Thus, in the r — > limit 
we must find also a unique prescription how to integrate lfM|l . 

The same reasoning can be carried out for individual particle trajectories and particle pair sep- 
aration, and since we have mostly developed formulae for the latter, we will use that setting. From 
the development of r for a small time-step At, as obtained in subsection IH.4L we could deduce that 
in the r — > limit our formulae lead to the Ito interpretation of (|Hn|). 

We can also check on the simulations (done for small enough time-step that we don't even need 
the sophistications of subsection 16. 4p that this is the case. Our simulation of particle separation is 
equivalent to tracing a single particle in a linear (in the spatial variable) velocity field, that is C/ = aR 
recalling Q. Hence the Lyapunov exponent of pair separation of the passive tracer is described by 
lf6?Hl with the identification R = f and U = v. Now the Lyapunov exponent obtained from ifM]) 
depends on the convention we adopt for the latter. Taking the Ito convention leads to substituting 
dr = dSr (recall (jHl) into l(5T1) and we get 

/ d log \ 1 1 

( ) = -^^CikjirirkVjri + -^Cik^urtri = Di{d - Ap){d - 1) 

In the simulations this is the value we get in the r — > limit. 

In comparison, Stratonovich interpretation of H6()p would require considering dr^ = dSijVj + 
{l/2){dSikdSki)ri = dSijrj + {l/2)Cik,kiri = dSijrj+Dip{d+2){d-l)ri which would give a Lyapunov 
exponent shifted by the non-zero Dip{d+ 2){d — 1). 

The fact that for inertial particles the Lagrangian equation (|6?Hl must always be taken within 
the Ito convention is due to the fundamental time-irreversible behavior of such particles. Indeed, 
particles have some knowledge of fiuid velocity in the past, which they forget only after about a time 
of the order of the Stokes time r, however they know nothing of the future. It is this smoothing 
of the fluid velocity at time-scale r, using only past events, that intuitively explains why the r — > 
limit leads invariably to a passive tracer which obeys ifM]) with Ito convention. 
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C Stratonovich convention for linearized velocity? 



This is just a short note to explain why one cannot naively first linearize (|6?Hl (in the position variable) 
and then change reading convention, say from Ito to Stratonovich. The two operations just do not 
commute with one another, and the physically meaningful order is the other one: linearizing the 
effective evolution of particles. 

To put very simply what fails, note that when passing from Ito to Stratonovich interpretation in 
(EOl, an additional drift term of the form {U .V)U appears: the Stratonovich interpretation of ((?)77|) 
corresponds to the Ito SDE 

AR{t) = dUt{R{t)) + ]^{{dUt.V)dUt) {R{t)) (61) 

Clearly, linearization of this term involves linearization of the derivative of dC/j, that comes from 
the quadratic part of dC/t, which would have been lost if we were to linearize dUt first and change 
convention after. 

Coming back to the simpler notation U , we can write the gradient of {U.V)U, in coordinate 
notation: 

dk[{Ujdj)Ui] = {dkUj){djUi) + UjdjdkUi = aijOjk + UjdkOij 

If the velocity field U has spatially homogeneous statistics, then {{Ujdj)Ui) is independent of position 
so its gradient is 0, implying that the Stratonovich term vanishes when we linearize JEB; we also 
have the identity: 

{(Tijcrjk) + {UjdkCFij) = 

Had we done linearization of the velocity field first, that is supposing a to be constant, we would 
have dkfJij = 0, and the linearization of the Stratonovich term would, incorrectly, reduce to the 
non-zero quantity \rk{(TijCFjk) proportional to (with contraction on repeated indices) ^Cij^k^k = 
2Di{d + 2){d — l)pri, vanishing only in the incompressible case p = 0. 

D Positive even order moments of pair separation 
D.l Arbitrary dimensions 

Let us rewrite equation (|T8|l as two coupled first-order equations for ip and ip = dip/dt in coordinate 
notation: 

dipi = ipidt , dtpi = -Eipidt + T~^{dSij)iij 

(recall definition of E from (|T7|) and of dS from ©) . The noteworthy fact is that this is a homogeneous 
equation in -0, so that moments of the same positive total order in the components will verify a 
closed system of equations: 

i,j=l k i^k j 

+E^^'^^r~'(d^/c)nv'rn^?) 

k i jT^k 

+E^^^"%^^r"'(d4)(d4)nc'n^?) 

k i jj^k 

+E^™^^'^r"Vr'"'(dV'fc)(dV'onc' n ^r) 
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which, upon replacing the dip and the dtp by their expressions, leads to 

k i^k j 

- E ^nfc(^f - n n i"') 

fc^!^; p,5 i j^k,l 

Introduce the multi-coefficients m and n for the mi and n^, and use the notation for adding 
1 to the k^^ coefficient of m, and for substracting 1 from the k*^ coefficient of m, etc. Also 
denote Cm,n = (Ilf j=i V'l^'V'j^)- Then we have the closed system for the Crn,n with [m] + [n] constant 
(brackets stand for the size of the multi-coefficient): 

k k 

+ E E "''^"^''""^V '^Cfcp^fcgCmP^g.nfc fc + E E '^kniT^'^Ckp,lqCrnP'i,n^ i (62) 

k p,q k^l p,q 

The idea which permits to compute the growth rate of ji/'p^ for G N is that it should be given 
by the topmost eigenvalue (Lyapunov exponent) of considered as a system of linear differential 
equations on Crn,n for ^nii + "^rii = 2N. The system in question is linear (in the Cm,n) with 
constant coefficients (i.e. independent of t), so it has an associated matrix, and what we need is the 
top eigenvalue of this matrix. In general this eigenvalue should be the exponential growth rate of 




D.2 Simplified version in 2D 

In the two-dimensional case a somewhat simpler formulation of the above may be given. Let us 
introduce tp = e*/^'^r, ip = dip/dt and their conjugates tp and ip. Recall 

^ +Vip = EtP 



Furthermore, denoting 

the equivalent of the system (|62)l is 
'^^ki 



dt2 



- kcl_^^i - E{n - kyi^^^i + lcli_^ - E{n - 1)4^1^^ 
+ {n - k){n - l){f3L + Pn)4+i,i+i 



This system has been derived in [TT] for similar purposes as here, in the framework of the Anderson 
localization problem. 



28 



References 



[1] Milton Abramowitz and Irene A. Stegun, editors. Handbook of mathematical functions with 
formulas, graphs, and mathematical tables. Dover Publications Inc., New York, 1992. Reprint 
of the 1972 edition. 

[2] J. Bee, A. Celani, M. Cencini, and S. Musacchio. Clustering and collisions of heavy particles in 
random smooth flows. Preprint at arXiv:iilin. CD/0407013. 

[3] G. Falkovich, A. Fouxon, and M. G. Stepanov. Acceleration of rain initiation by cloud turbulence. 
Nature, 419:151-4, 2002. 

[4] G. Falkovich, K. Gaw§dzki, and M. Vergassola. Particles and fields in fluid turbulence. Rev. 
Modern Phys., 73(4):913-975, 2001. 

[5] A. Fouxon and P. Horvai. Inertial particles and the Anderson localization problem. In consid- 
eration for Phys. Rev. Lett., 2006. Preprint will be in arXiv: cond-mat. 

[6] Bertrand I. Halperin. Green's functions for a particle in a one-dimensional random potential. 
Phys. Rev. (2), 139:A104-A117, 1965. 

[7] M. R. Maxey and Riley J. Equation of motion of a small rigid sphere in a nonuniform flow. 
Phys. Fluids, 26:883-889, 1983. 

[8] B. Mehlig, M. Wilkinson, K. Duncan, T. Weber, and Ljunggren M. On the aggregation 
of inertial particles in random flows. Submitted to Phys. Rev. E (3), 2005. Preprint at 
arXiv : cond-inat/0410518. 

[9] Bernhard Mehlig and Michael Wilkinson. Coagualtion by random velocity fields as a Kramers 
problem. Phys. Rev. Lett., 92(25) :250602-l-4, 2004. 

[10] Leonid I. Piterbarg. The top Lyapunov exponent for a stochastic flow modeling the upper ocean 
turbulence. SI AM J. Appl. Math., 62(3):777-800 (electronic), 2001/02. 

[11] H. Schomerus and M. Titov. Statistics of finite-time Lyapunov exponents in a random time- 
dependent potential. Phys. Rev. E (3), 66(6):066207, 11, 2002. 



29 




Figure 2: Adimensionalized Lyapunov exponent in function of Stokes number for different values of p. 
We see perfect agreement between theory (solid lines only, i.e. those with the empty boxes) and numerics 
(boxes, solid or empty) for <p = —1/2 and p = 3/2. 
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Figure 3: Relative error of numerical results with respect to the formula conjectured in [10] and reproduced 
in our paper as formula (f57jl . for p < 1/2. A zoom (not represented) on the curves for very small St is 
not incompatible with the prediction (according to the divergent asymptotic expansion in [H]) that all 
derivatives of the relative error curves should vanish at St = 0, but quality of our data in that range is 
too poor. 





0.5 









-0.5 








-1 


(N 






-1.5 


T-H 




-2 




-2.5 




-3 



1 \ r 







0' 



0' 







\ \ L 



1 \ r 








1 r 
-e 



J \ L 



-<3> 



p=\ --e- 

^ = 3/2 --<3>- 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 



Figure 4: Non-collapse of numeric curves for p > 1/2, contradicting (f59ll 
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